#include "AnaDef.h"
#include "XGLUtils.h"
#include "style.h"

const Double_t EPSILON = 1e-12;

const Int_t fcol[3]={kRed,  kBlue, kGreen+3};
const Int_t fmsty[3]={27,24,30};//20,21,24};//,25 };
const Int_t flsty[3]={kSolid, kDotted,   kDashed};
const Int_t ffsty[3]={0,0,0};//1001,1001,1001};

Int_t fgrid = -999;
Int_t fmode = -999;

const Int_t fnjpt = 3;
const TString sjpt[]={"5-10", "10-15", "15-20"};

const Int_t fntype = 3;
const TString stype[]={"p+#bar{p}","#pi^{+}+#pi^{-}","K^{+}+K^{-}"};

const Int_t fntheo = 3;
const TString stheo[]={"Perugia0","Perugia2010NoCR","Perugia2011"};

TList *flist = new TList;

const Int_t fillcol[]={kGray+2,kGray+2, kGray+2};

//===========================================================

const Int_t nGr(){return fnjpt*fntype;}

Double_t grMax(const TGraphAsymmErrors *gr)
{
  Double_t gm = -1e10;
  for(Int_t ii=0; ii<gr->GetN(); ii++){
    if(gr->GetY()[ii]>gm){
      gm = gr->GetY()[ii];
    }
  }
  return gm;
}

Int_t IDx()
{
  const Int_t idx = fgrid/10;
  if(idx!=1 && idx!=2){
    printf("bad idx %d %d\n", fgrid, idx); exit(1);
  }
  return idx;
}

Int_t IDy()
{
  const Int_t idy= fgrid%10;
  if(idy!=6 && idy!=7){
    printf("bad idy %d %d\n", fgrid, idy); exit(1);
  }
  return idy;
}


Int_t getJpt(TString gn)
{
  gn.ReplaceAll("_","-");

  Int_t ijpt = -999;
  for(Int_t ii = 0; ii<fnjpt; ii++){
    if(gn.Contains(sjpt[ii])){
      if(ijpt>=0){
        printf("getJpt bad gn %s %d\n", gn.Data(), ijpt); exit(1);
      }
      ijpt = ii;
    }
  }

  if(ijpt<0){
    printf("getJpt bad again! gn %s %d\n", gn.Data(), ijpt); exit(1);
  }

  return ijpt;
}


Int_t getType(const TString gn)
{
  Int_t itype = -999;
  for(Int_t ii=0; ii<fntype; ii++){
    if(gn.Contains(Form("grstat%d",ii)) || gn.Contains(Form("grsys%d",ii))){

      if(itype>=0){
        printf("getType bad gn %s %d\n", gn.Data(), itype); exit(1);
      }

      itype=ii;
    }
  }
  if(itype<0){
    printf("getType bad again gn %s %d\n", gn.Data(), itype); exit(1);
  }
  return itype;
}


Int_t getTheo(TString gn)
{
  gn.ReplaceAll("-","");
  gn.ReplaceAll("_","");

  Int_t itheo=-999;
  for(Int_t ii=0; ii<fntheo; ii++){
    if(gn.Contains(stheo[ii])){

      if(itheo>=0){
        printf("getTheo bad gn %s %d\n", gn.Data(), itheo); exit(1);
      }

      itheo=ii;
    }
  }

  if(itheo<0){
    //can happen for read data
    //printf("getTheo bad again! gn %s %d\n", gn.Data(), itheo); exit(1);
  }

  return itheo;
}

void safePlot(TGraphAsymmErrors* gr, const TString opt)
{
  style::ResetStyle(gr);

  const TString gn=gr->GetName();

  const Int_t itype = getType(gn);
  const Int_t ijpt = getJpt(gn);
  const Int_t itheo = getTheo(gn);
  const TString ytit=gr->GetYaxis()->GetTitle();

  //===================================================

  Int_t id = ijpt;
  if(itheo>=0)
    id = itheo;

  gr->SetLineColor(fcol[id]);
  gr->SetMarkerColor(fcol[id]);
  gr->SetMarkerStyle(fmsty[id]);
  if(itheo>=0){
    gr->SetLineStyle(flsty[itheo]);
  }

  if(gn.Contains("sys")){
    gr->SetFillColor(fillcol[id]);
    gr->SetFillStyle(ffsty[id]);
    gr->SetLineColor(fcol[id]);
  }
  
  if(fmode==3 || ytit.Contains("MC/Data")){
    if(itheo<0){
      gr->SetFillStyle(ffsty[id]);
      gr->SetFillColor(fillcol[id]);
      gr->SetMarkerColor(fillcol[id]);
      gr->SetMarkerStyle(1);
      gr->SetMarkerSize(0);
      gr->SetLineColor(fillcol[id]);
    }
    else{
      if(gr->GetLineStyle()!=kSolid){
        gr->SetLineWidth(3);
      }
    }
  }
  
  //==================================================
  if(IDx()==1){
    if(IDy()==6){
      gr->GetXaxis()->SetLimits(0.12, 25);
    }
    else{
      gr->GetXaxis()->SetLimits(0.12, 25);//lower limit > 0.2 too crowded
    }
  }
  else{
    gr->GetXaxis()->SetLimits(-0.05,1.05);
  }

  gr->GetXaxis()->SetMoreLogLabels();
  gr->GetYaxis()->SetNdivisions(505);

  //====================================================
 
  if(IDx()==1){
    gr->GetXaxis()->SetTitle("#it{p}_{T}^{track} (GeV/#it{c})");
  }
  else if(IDx()==2){
    gr->GetXaxis()->SetTitle("#it{z}^{ch}");
  }

  if(itype!=0){
    gr->GetXaxis()->SetTitle("");
  }

  if(fmode!=1){
  //also for maxtrix comparison
    if(ijpt!=0){
      gr->GetYaxis()->SetTitle("");
    }
  }
  //==================================================
  
  flist->Add(gr);

  XGLUtils::PrintGr("printbeforedraw", gr);

  gr->Draw(opt); printf("printatdraw: %s color %d %d marker %d title ploted %s title old %s\n", gr->GetName(), gr->GetLineColor(), gr->GetMarkerColor(), gr->GetMarkerStyle(), gr->GetYaxis()->GetTitle(), ytit.Data());

}

Int_t grIter(const Int_t itype, const Int_t ijpt)
{
  return itype*fnjpt + ijpt;
}

Int_t theoIter(const Int_t itheo, const Int_t igr)
{
  return itheo*nGr() + igr;
}

void ReadData(TGraphAsymmErrors * grstat[], TGraphAsymmErrors * grsys[])
{
  for(Int_t ijpt=0; ijpt<fnjpt; ijpt++){
    TFile * tmpf=new TFile(Form("data/JPT%s_GrID%d.root", sjpt[ijpt].Data(), fgrid));
    if(!tmpf->IsOpen()){
      printf("file not open %s\n", tmpf->GetName()); exit(1);
    }

    TString nsj(sjpt[ijpt]);
    nsj.ReplaceAll("-","_");

    TString stmp;
    for(Int_t itype=0; itype<fntype; itype++){
      stmp=Form("grstat%d", itype);
      const Int_t igr = grIter(itype, ijpt);
      grstat[igr] = (TGraphAsymmErrors *) tmpf->Get(stmp);
      if(!grstat[igr]){
        printf("no grstat itype %d ijpt %d igr %d %s\n", itype, ijpt, igr, stmp.Data());
        tmpf->ls();
        exit(1);
      }
      grstat[igr]->SetName(Form("data_%s_%s", nsj.Data(), grstat[igr]->GetName()));

      if(IDy()==6){
        if(IDx()==1){
          grstat[igr]->GetYaxis()->SetTitle(Form("1/N_{jets}dN/d#it{p}_{T}^{track}(GeV/#it{c})^{-1}"));
        }
        else if(IDx()==2){
          grstat[igr]->GetYaxis()->SetTitle(Form("1/N_{jets} dN/d#it{z}^{ch}"));
        }
      }
      else{
        grstat[igr]->GetYaxis()->SetTitle("particle ratio");
      }
      //===========================
      
      stmp=Form("grsys%d", itype);
      grsys[igr] = (TGraphAsymmErrors *) tmpf->Get(stmp);
      if(!grsys[igr]){
        printf("no grsys itype %d ijpt %d igr %d %s\n", itype, ijpt, igr, stmp.Data());
        tmpf->ls();
        exit(1);
      }
      grsys[igr]->SetName(Form("data_%s_%s", nsj.Data(), grsys[igr]->GetName()));
      grsys[igr]->GetYaxis()->SetTitle(grstat[igr]->GetYaxis()->GetTitle());
    }
  }
}

void ReadTheo(TGraphAsymmErrors *grtheo[])
{
  for(Int_t itheo=0; itheo<fntheo; itheo++){
    for(Int_t ijpt=0; ijpt<fnjpt; ijpt++){
      TFile * tmpf=new TFile(Form("theory/outTheoryCombined_%s_JPT%s_GrID%d.root", stheo[itheo].Data(), sjpt[ijpt].Data(), fgrid));
      if(!tmpf->IsOpen()){
        printf("file not open %s\n", tmpf->GetName()); exit(1);
      }

      TString stmp;
      for(Int_t itype=0; itype<fntype; itype++){
        stmp=Form("grstat%d", itype);
        const Int_t igr = theoIter(itheo, grIter(itype, ijpt));
        grtheo[igr] = (TGraphAsymmErrors *) tmpf->Get(stmp);
        if(!grtheo[igr]){
          printf("no grstat itheo %d itype %d ijpt %d igr %d %s\n",itheo, itype, ijpt, igr, stmp.Data());
          tmpf->ls();
          exit(1);
        }
        TString nsj(sjpt[ijpt]);
        nsj.ReplaceAll("-","_");
        grtheo[igr]->SetName(Form("%s_%s_%s", stheo[itheo].Data(), nsj.Data(), grtheo[igr]->GetName()));
      }
    }
  }

}

void clearTheo(TGraphAsymmErrors * grs[])
{
  Double_t xcutoff = -999;
  for(Int_t itheo=0; itheo<fntheo; itheo++){
  for(Int_t itype=0; itype<fntype; itype++){
    if(IDx()==1){
      if(itype==1){
        xcutoff = 0.15;
      }
      else{
        xcutoff = 0.35;
      }
    }

    for(Int_t ijpt = 0; ijpt<fnjpt; ijpt++){
      const Int_t ii = theoIter(itheo, grIter(itype, ijpt));
      
      TGraphAsymmErrors* gr = grs[ii];
      
      //----------------                                                                                                                            
      const Double_t oldlastx0 = gr->GetX()[gr->GetN()-1];
      
      Int_t istart = -999;
      for(Int_t ii=0; ii<gr->GetN(); ii++){
        if(gr->GetX()[ii]>xcutoff){
          istart = ii;
          break;
        }
      }
      if(istart<0){
        printf("clearTheo istart<0\n");
        for(Int_t ii = 0; ii < gr->GetN(); ii++){
          printf("raw %d/%d x %.12e %e\n", ii, gr->GetN(), gr->GetX()[ii], gr->GetY()[ii]);
        }
        exit(1);
      }
      
      const Int_t nn = gr->GetN()-istart;
      for(Int_t ii=0; ii<nn; ii++){
        const Int_t iold = ii+istart;
        const Double_t xx = gr->GetX()[iold];
        const Double_t yy = gr->GetY()[iold];
        const Double_t exl = gr->GetEXlow()[iold];
        const Double_t exh = gr->GetEXhigh()[iold];
        const Double_t eyl = gr->GetEYlow()[iold];
        const Double_t eyh = gr->GetEYhigh()[iold];
        
        gr->SetPoint(ii, xx, yy);
        gr->SetPointError(ii, exl, exh, eyl, eyh);
      }
      gr->Set(nn);
      
      const Double_t newlastx0 = gr->GetX()[gr->GetN()-1];
      if(fabs(oldlastx0-newlastx0)>EPSILON){
        printf("clearTheo old != new %e %e\n", oldlastx0, newlastx0); exit(1);
      }
      //-------------------------                                                                                                                     
      printf("checkfirstpoint %s %e\n", gr->GetName(), gr->GetX()[0]);
    }
  }
  }
}

TGraphAsymmErrors *scaleGr(const TGraphAsymmErrors *gtar, const TGraphAsymmErrors *gref)
{
  //Int_t style::GetCommonLimits(const Int_t ngr, const TGraphAsymmErrors * grs[], Int_t id0s[], Int_t id1s[])

  const Int_t ngr = 2;
  Int_t id0s[ngr];
  Int_t id1s[ngr];
  const TGraphAsymmErrors * grs[]={gtar, gref};
  const Int_t npt = style::GetCommonLimits(ngr, grs, id0s, id1s);
  if(npt<10){
    printf("scaleGr wrong npt %d\n", npt);exit(1);
  }

  //copy number of points and xarray from gREF, with both name
  const TString nref = gref->GetName();
  const TString ntar = gtar->GetName();
  TGraphAsymmErrors * gout = new TGraphAsymmErrors(npt);
  gout->SetName(Form("%sscaled%s", nref==ntar?"self":"", gtar->GetName()));

  for(Int_t iout=0; iout<npt; iout++){
    const Int_t iref=id0s[1]+iout;
    const Int_t itar=id0s[0]+iout;
    
    const Double_t xx = gref->GetX()[iref];

    if(fabs(xx-gtar->GetX()[itar])>EPSILON){
      printf("bad point ref %d %e tar %d %e\n", iref, xx, itar, gtar->GetX()[itar]); exit(1);
    }
    
    const Double_t yref = gref->GetY()[iref];
    gout->SetPoint(iout, xx,   gtar->GetY()[itar]/yref);
    gout->SetPointError(iout, gtar->GetEXlow()[itar], gtar->GetEXhigh()[itar], gtar->GetEYlow()[itar]/yref, gtar->GetEYhigh()[itar]/yref);
  }

  return gout;
}

void cookRatData(TGraphAsymmErrors *grRat[], TGraphAsymmErrors *grold[])
{
  for(Int_t itype=0; itype<fntype; itype++){
    for(Int_t ijpt=0; ijpt<fnjpt; ijpt++){
      const Int_t igr=grIter(itype, ijpt);
      const Int_t jr = fgrid==16?2:0;
      const Int_t iref = grIter(itype, jr);
      grRat[igr]=scaleGr(grold[igr], grold[iref]);
      grRat[igr]->GetYaxis()->SetTitle(Form("#splitline{scaled to}{%s GeV/#it{c}}", sjpt[jr].Data()));
    }
  }
}

void cookScData(TGraphAsymmErrors *grSc[], TGraphAsymmErrors *grold[])
{
  for(Int_t igr=0; igr<nGr(); igr++){
    grSc[igr]=scaleGr(grold[igr], grold[igr]);
    //grSc[igr]->GetYaxis()->SetTitle(Form("%s MC/Data", grold[igr]->GetYaxis()->GetTitle()));
    grSc[igr]->GetYaxis()->SetTitle("MC/Data");
  }
}

void cookScTheo(TGraphAsymmErrors *grScTheo[], TGraphAsymmErrors *grtheo[], TGraphAsymmErrors *grstat[])
{
  for(Int_t itheo=0; itheo<fntheo; itheo++){
    for(Int_t igr=0; igr<nGr(); igr++){
      const Int_t id = theoIter(itheo, igr);
      grScTheo[id] = scaleGr(grtheo[id], grstat[igr]);
      //grScTheo[igr]->GetYaxis()->SetTitle(Form("%s MC/Data", grstat[igr]->GetYaxis()->GetTitle()));
      grScTheo[id]->GetYaxis()->SetTitle("MC/Data");
    }
  }
}

void PlotJetEvol(const Bool_t kzoom, TGraphAsymmErrors *grstat[], TGraphAsymmErrors*grsys[], const Int_t itype)
{
  const Int_t jptid[]={2,1,0};
  for(Int_t rawi=0; rawi<fnjpt; rawi++){
    const Int_t ijpt = IDx()==1?jptid[rawi]:rawi;
    const Int_t gid = grIter(itype, ijpt);

    const Double_t zmin = -0.4;
    const Double_t zmax = 2.4;
    if(IDy()==6){
      if(IDx()==1){
        if(!kzoom){
          grstat[gid]->SetMinimum(3e-6);
          grstat[gid]->SetMaximum(8);
        }
        else{
          grstat[gid]->SetMinimum(zmin);
          grstat[gid]->SetMaximum(zmax);
        }
      }
      else{
        if(!kzoom){
          grstat[gid]->SetMinimum(4e-3);
          grstat[gid]->SetMaximum(50);
        }
        else{
          grstat[gid]->SetMinimum(zmin);
          grstat[gid]->SetMaximum(zmax);
        }
      }
    }
    else{
      if(IDx()==1){
        if(!kzoom){
          grstat[gid]->SetMinimum(-0.03);
          grstat[gid]->SetMaximum(0.8);
        }
        else{
          grstat[gid]->SetMinimum(0.25);
          grstat[gid]->SetMaximum(zmax);
        }
      }
      else{
        if(!kzoom){
          grstat[gid]->SetMinimum(-0.02);
          grstat[gid]->SetMaximum(0.7);
        }
        else{
          grstat[gid]->SetMinimum(0.25);
          grstat[gid]->SetMaximum(zmax);
        }
      }
    }
    grsys[gid]->SetMinimum(grstat[gid]->GetMinimum());
    grsys[gid]->SetMaximum(grstat[gid]->GetMaximum());
    if(!kzoom){    
      safePlot(grsys[gid],rawi?"2":"a2");     
    }
    else{
      //safePlot(grsys[gid],rawi?"3":"a3");     
      safePlot(grsys[gid],rawi?"2":"a2");     
    }
    safePlot(grstat[gid], rawi?"p":"p"); 

  }

}

void PlotTheo(const Bool_t kratio, TGraphAsymmErrors *grScStat[], TGraphAsymmErrors *grScSys[], TGraphAsymmErrors *grScTheo[], const Int_t itype, const Int_t ijpt)
{
  const Int_t igr = grIter(itype, ijpt);
  Double_t hmin = 0.3;
  Double_t hmax = 2.8;
  if(!kratio){
    if(IDy()==6){
      if(gPad->GetLogy()){
        if(IDx()==1){
          hmin = 3e-6;
          hmax = 8;
        }
        else{
          hmin = 4e-3;
          hmax = 50;
        }
      }
      else{//just see proton at 5-10
        hmin = -0.05*grMax(grScStat[0]);
        hmax = grMax(grScStat[0])*1.5;
      }
    }
    else{
      if(IDx()==1){
        hmin = -0.03;
        hmax = 0.77;
      }
      else{
        hmin = -0.03;
        hmax = 0.77;
      }
    }
  }
  
  grScStat[igr]->SetMinimum(hmin);
  grScStat[igr]->SetMaximum(hmax);
  grScSys[igr]->SetMinimum(hmin);
  grScSys[igr]->SetMaximum(hmax);
  safePlot(grScSys[igr],"a2");
  //safePlot(grScSys[igr],"a3");
  safePlot(grScStat[igr],"p");

  const Double_t lastx = grScStat[igr]->GetX()[grScStat[igr]->GetN()-1];

  for(Int_t itheo=0; itheo<fntheo; itheo++){
    const Int_t thit = theoIter(itheo, igr);
    grScTheo[thit]->SetMinimum(hmin);
    grScTheo[thit]->SetMaximum(hmax);

    //don't draw overflow
    Int_t ilast = -999;
    for(Int_t ipt = 0; ipt<grScTheo[thit]->GetN(); ipt++){
      if( fabs(grScTheo[thit]->GetX()[ipt]-lastx)<EPSILON ){
        ilast = ipt;
        break;
      }
    }
    if(ilast<0){
      printf("plotTheo bad no ilast\n"); exit(1);
    }
    grScTheo[thit]->Set(ilast+1);

    safePlot(grScTheo[thit], itheo?"lX":"lX");
  }
}

void DrawParticle(const Int_t ipad, const Int_t nx, const Int_t itype, const Double_t yscale=1,  Double_t tx = -999, Double_t ty = -999 )
{
  if(ipad<nx)
    return;

  if(tx<-990 && ty<-990){
    tx = 0.78;
    ty = 0.85;
    if(IDy()==7 &&  yscale>0.99){
      ty=0.8;
      if(ipad==nx){
        if(IDx()==1){
          ty=0.15;
        }
        else{
          tx=0.25;
        }
      }
    }
  }
  
  TLatex *lpar=new TLatex(tx, ty*yscale, IDy()==6?stype[itype]:Form("#frac{%s}{%s}",stype[itype].Data(),stype[1].Data()));
  style::ResetStyle(lpar);
  lpar->Draw();
}


void DrawJet(const Int_t ipad, const Int_t nx, const Int_t ijpt, Double_t tx = -999, Double_t ty = -999, Double_t ltscale=-999)
{
  if(ipad%nx!=0)
    return;

  if(tx<-990)
    tx = 0.25;

  if(ty<-990)
    ty = 0.8;
  
  TH1D *hdum=new TH1D;
  hdum->SetFillStyle(ffsty[2]);
  hdum->SetFillColor(fillcol[2]);
  hdum->SetLineColor(kBlack);
  hdum->SetMarkerColor(kBlack);
  hdum->SetMarkerStyle(25);
  hdum->SetMarkerSize(1);

  //printf("testtestty %e %e\n", tx, ty);

  //TLegend * lg=new TLegend(tx,ty,  tx+0.25, ty+0.15);
  //lg->AddEntry(hdum, sjpt[ijpt]+" GeV/#it{c}","lf");
  //lg->SetHeader(sjpt[ijpt]+" GeV/#it{c}");
  //style::ResetStyle(lg);
  //lg->SetBorderSize(1);
  //  lg->Draw();

  //TLatex *lpar=new TLatex(tx, ty, sjpt[ijpt]+" GeV/#it{c}");
  TLatex *lpar=new TLatex(tx, ty, Form("#it{p}_{T,jet}^{ch} %s GeV/#it{c}", sjpt[ijpt].Data()));

  if(ltscale<-990){
    ltscale = 0.8;
  }

  style::ResetStyle(lpar,ltscale);
  lpar->Draw();
}

void DrawAna(const Int_t ipad, const Int_t nx)
{
  if( ipad!=nx+1 )
    return;

  //const TString jeth0("#splitline{FastJet anti-#it{k}_{T}}{#it{R}=0.4; |#eta^{jet}|<0.5}");
  const TString jeth0("#splitline{anti-#it{k}_{T}}{#it{R}=0.4; |#eta^{jet}|<0.5}");
  const TString jeth1("#splitline{#it{p}_{T}^{track}>0.15 GeV/#it{c}}{|#eta^{track}|<0.9}");
  
  TLatex *lt0=0x0;
  TLatex *lt1=0x0;
  
  const Double_t x0 = 0.05;

  if(IDy()==7){
    lt0=new TLatex(x0,0.82, jeth0); 
    lt1=new TLatex(x0,0.62, jeth1);
  }
  else{
    lt0=new TLatex(x0,0.3, jeth0); 
    lt1=new TLatex(x0,0.1, jeth1);
  }
  
  style::ResetStyle(lt0,0.8); 
  style::ResetStyle(lt1,0.8);
  lt0->Draw();
  lt1->Draw();
}


void Draw2Ana(const Int_t ipad, const Int_t nx, Double_t x0=-999, Double_t y0 = -999)
{
  if( ipad!=nx+1 )
    return;

  if(y0<-990)
    y0=0;

  //const TString jeth0("FastJet anti-#it{k}_{T}; #it{R}=0.4; |#eta^{jet}|<0.5");
  const TString jeth0("anti-#it{k}_{T}; #it{R}=0.4; |#eta^{jet}|<0.5");
  const TString jeth1("#it{p}_{T}^{track}>0.15 GeV/#it{c}; |#eta^{track}|<0.9");
  
  TLatex *lt0=0x0;
  TLatex *lt1=0x0;
  
  if(x0<-990)
    x0 = IDy()==6?0.05:0.35;

  const Double_t yscale = IDy()==6? 0.02 : 0;
  lt0=new TLatex(x0,y0+0.88-yscale, jeth0); 
  lt1=new TLatex(x0,y0+0.78-yscale, jeth1);
 
  style::ResetStyle(lt0,0.7); 
  style::ResetStyle(lt1,0.7);
  lt0->Draw();
  lt1->Draw();
}

void DrawPre(const Int_t ipad, const Int_t nx)
{
 if( ipad!=nx )
   return;

 TLegend * lg = 0x0;
 if(IDy()==6){
   lg = new TLegend(0.25, 0.02, 0.5, 0.40);
 }
 else{
   if(IDx()==1){
     lg = new TLegend(0.25, 0.52, 0.5, 0.90);
   }
   else{
     lg = new TLegend(0.6, 0.02, 0.85, 0.40);
   }
 }

 for(Int_t dummy=0; dummy<fnjpt; dummy++){
   TH1D *hdum=new TH1D;
   hdum->SetLineColor(fcol[dummy]);
   hdum->SetMarkerStyle(fmsty[dummy]);
   hdum->SetMarkerColor(fcol[dummy]);
   hdum->SetFillColor(fcol[dummy]);
   hdum->SetFillStyle(ffsty[dummy]);
   lg->AddEntry(hdum,sjpt[dummy]+" GeV/#it{c}", "lp");//"lfp");
   //lg->AddEntry(hdum,Form("p_{T,jet}^{ch} %s GeV/#it{c}", sjpt[dummy].Data()), "lfp");
 }
 lg->SetHeader("#splitline{pp #sqrt{s} = 7 TeV}{ALICE Preliminary}");
 style::ResetStyle(lg,0.2);
 lg->Draw();
}

void Draw2Pre(const Int_t ipad, const Int_t nx, Double_t x0 = -999, Double_t y0 = -999, Double_t y2=-999)
{
  if(x0<-990)
    x0=IDy()==6?0.05:0.35;

  if(y0<-990)
    y0 = 0.5;

  if( ipad==nx+1 ){
    TLegend * lg = 0x0;
    lg = new TLegend(x0, y0, 0.5, y0+0.4);

    TH1D *hdum=new TH1D;
    hdum->SetFillStyle(ffsty[2]);
    hdum->SetFillColor(fillcol[2]);
    hdum->SetLineColor(kBlack);
    hdum->SetMarkerColor(kBlack);
    hdum->SetMarkerStyle(25);
    hdum->SetMarkerSize(1);

    lg->AddEntry(hdum, "data","lf");

    for(Int_t dummy=0; dummy<fnjpt; dummy++){
      TH1D *hdum=new TH1D;
      hdum->SetLineColor(fcol[dummy]);
      hdum->SetLineStyle(flsty[dummy]);
      if(flsty[dummy]!=kSolid){
        hdum->SetLineWidth(2);
      }
      hdum->SetMarkerStyle(fmsty[dummy]);
      hdum->SetMarkerColor(fcol[dummy]);
      TString tmpent(stheo[dummy]);
      tmpent.ReplaceAll("2010NoCR","0NoCR");
      lg->AddEntry(hdum,Form("PYTHIA %s", tmpent.Data()),"l");
    }
    style::ResetStyle(lg,0.4);
    lg->Draw();
  }
  else if(ipad==2*nx+1){
     TLegend * lg = 0x0;
     if(y2>-990)
       y0=y2;

     lg = new TLegend(x0, y0+0.1, 0.5,  y0+0.5);

     lg->SetHeader("#splitline{pp #sqrt{s} = 7 TeV}{ALICE Preliminary}");
     style::ResetStyle(lg);
     lg->Draw();
  }
}

int main(int argc, char * argv[])
{
  for(Int_t ii=0; ii<argc; ii++){
    printf("%d: %s\n", ii, argv[ii]);
  }

  Bool_t kflipLogy = kFALSE;
  if(argc==5){

  }
  else if(argc==4){
    kflipLogy = kTRUE;
  }


  if(argc!=3){
    printf("argc!=3\n"); return 1;
  }

  fgrid = atoi(argv[1]);
  fmode=atoi(argv[2]);

  //_______________________________________________
  TGraphAsymmErrors * grstat[nGr()];
  TGraphAsymmErrors * grsys[ nGr()];
  ReadData(grstat, grsys);

  //_______________________________________________
  TGraphAsymmErrors * grRatStat[nGr()];
  TGraphAsymmErrors * grRatSys[ nGr()];
  cookRatData(grRatStat, grstat);
  cookRatData(grRatSys, grsys);

  //_______________________________________________
  TGraphAsymmErrors * grScStat[nGr()];
  TGraphAsymmErrors * grScSys[nGr()];
  cookScData(grScStat, grstat);
  cookScData(grScSys, grsys);
  //_______________________________________________
  TGraphAsymmErrors * grtheo[fntheo*nGr()];
  ReadTheo(grtheo);

  //_______________________________________________
  TGraphAsymmErrors * grScTheo[fntheo*nGr()];
  cookScTheo(grScTheo, grtheo, grstat);

  //_______________________________________________
  //_______________________________________________
  style::SetGlobalStyle(0);

  //pion, kaon, proton
  Int_t *parid=0x0;
  if(IDy()==6){
    parid=new Int_t[3];
    parid[0]=1;
    parid[1]=2;
    parid[2]=0;
  }
  else{
    parid=new Int_t[2];
    parid[0]=2;
    parid[1]=0;
  }

  const Double_t cw = IDy()==6?1200:800;
  const Double_t ch = 450;
  const Int_t nx = IDy()==6? 3 : 2;

  Int_t npad = 0;

  style::fgkLineWidth=1;
  //_______________________________________________
  //_______________________________________________
  
  if(fmode==1){

  TCanvas * c1 = new TCanvas("IDFFpp7TeV_data", "", cw,ch);
  TList * list1 = new TList;
  npad = 0;
  style::DividePad(c1, nx, 4, 0.2, 0.01, 0.03, 0.4+(IDx()==1?0.05:0), list1, npad, 1);  

  for(Int_t ipad=0; ipad<npad; ipad++){
    TPad * pp = (TPad*)list1->At(ipad);
    pp->cd();

    if(IDx()==1){
      pp->SetLogx();
    }

    const Int_t itype = parid[ipad%nx];

    printf("ipad %d itype %d\n", ipad, itype);

    const Double_t rs = 0.15;

    if(ipad>=nx){
      const Double_t f3=0.55;
      style::fgkTitleSize = rs*f3;
      style::fgkTextSize = rs*f3;
      
      style::fgkYTitleOffset = 1.15;
      style::fgkTickLength = 0.02;

      //=============================
      if(IDy()==6){
        pp->SetLogy();
      }
      style::fgkMarkerSize= IDy()==7?1:0.5;
      PlotJetEvol(kFALSE, grstat, grsys, itype);
    }
    else{
      style::fgkTitleSize = rs;
      style::fgkTextSize = rs;
      style::fgkYTitleOffset = 0.6;
      style::fgkXTitleOffset = 1.3;
      style::fgkMarkerSize= 0.5;
      style::fgkTickLength = 0.02*2;

      PlotJetEvol(kTRUE, grRatStat, grRatSys, itype);
    }

    DrawPre(ipad,nx);
    DrawAna(ipad, nx);
    DrawParticle(ipad, nx, itype);
  }

  c1->Draw();
  c1->Print(Form("%s_ID%d.eps",c1->GetName(), fgrid));
  c1->Print(Form("%s_ID%d.root",c1->GetName(), fgrid));
  }
  else if(fmode==2){
  //_______________________________________________
  //_______________________________________________
 
  npad=0;
  TCanvas * c2 = new TCanvas("IDFFpp7TeV_MC", "", cw,ch);
  TList *list2 = new TList;
  style::DividePad(c2, nx, 3, 0.2, 0.01, 0.03, IDx()==2 ? 0.28: 0.4, list2, npad, 2);

  for(Int_t ipad=0; ipad<npad; ipad++){
    TPad * pp = (TPad*)list2->At(ipad);
    pp->cd();

    if(IDx()==1){
      pp->SetLogx();
    }

    const Int_t itype = parid[ipad%nx];

    printf("ipad %d itype %d\n", ipad, itype);

    const Double_t rs = 0.15;

    const Double_t ff = IDy()==7?0.8:1;
    if(ipad<nx){  
      style::fgkTitleSize = rs*0.8*ff;
      //style::fgkTextSize = rs*0.9*ff;
      style::fgkTextSize = rs*(IDx()==1?0.6:0.75)*ff;
      style::fgkYTitleOffset = 0.5/ff;
      style::fgkXTitleOffset = 1.3/ff*(IDx()==2?0.6:1);
    }
    else{
      style::fgkTitleSize = rs*1.*ff;
      style::fgkTextSize = rs*1.*ff;
      style::fgkYTitleOffset = 0.4/ff;
      style::fgkXTitleOffset = 1.3/ff;
    }
    
    const Int_t ijpt = 2-ipad/nx;
    printf("ipad %d ijpt %d\n", ipad, ijpt);
    if(ijpt<0||ijpt>2){
      printf("bad ijpt %d\n", ijpt);exit(1);
    }   
    
    PlotTheo(kTRUE, grScStat, grScSys, grScTheo, itype, ijpt);

    //==========================
    DrawParticle(ipad, nx*2, itype, 0.8);
    DrawJet(ipad,nx, ijpt, -999,-999, ipad<nx?0.8:0.8);
    Draw2Pre(ipad,nx);
    Draw2Ana(ipad,0);
  }

  c2->cd();
  c2->Draw();
  c2->Print(Form("%s_ID%d.eps",c2->GetName(), fgrid));
  c2->Print(Form("%s_ID%d.root",c2->GetName(), fgrid));
  }
  else if(fmode==3){
  //_______________________________________________
  //_______________________________________________
 
  npad=0;
  TCanvas * c3 = new TCanvas(Form("IDFFpp7TeV_comb%s",kflipLogy?"_WrongYScale":""), "", cw,ch*2);
  TList *list3 = new TList;
  style::DividePad(c3, nx, 3, /*IDy()==7? 0.15:*/ 0.25, 0.01, 0.03, IDx()==1?0.25:  0.2, list3, npad, 2);

  //clear
  if(IDx()==1){
    clearTheo(grtheo);
  }

  for(Int_t ipad=0; ipad<npad; ipad++){
    TPad * pp = (TPad*)list3->At(ipad);
    pp->cd();

    if(IDx()==1){
      pp->SetLogx();
    }

    if(IDy()==6){
      pp->SetLogy();
      if(kflipLogy){
        pp->SetLogy(0);
      }
    }

    const Int_t itype = parid[ipad%nx];

    printf("ipad %d itype %d\n", ipad, itype);

    const Double_t rs = 0.15;

    const Double_t ff = 0.5;
    if(ipad<nx){  
      style::fgkTitleSize = rs*1.*ff;
      style::fgkTextSize = rs*1*ff;
      style::fgkYTitleOffset = 0.5/ff*1.3;
      style::fgkXTitleOffset = 1.3/ff*0.4*(IDx()==1?1.5:1);
    }
    else{
      //style::fgkTitleSize = rs*1.5*ff;
      style::fgkTitleSize = rs*1.3*ff;
      style::fgkTextSize = rs*1.5*ff/1.2;
      //style::fgkYTitleOffset = 0.4/ff*1.3;
      style::fgkYTitleOffset = 0.4/ff*1.5;
      //style::fgkXTitleOffset = 1.3/ff;
      style::fgkXTitleOffset = 1.3/ff*1.2;
    }
    
    const Int_t ijpt = 2-ipad/nx;
    printf("ipad %d ijpt %d\n", ipad, ijpt);
    if(ijpt<0||ijpt>2){
      printf("bad ijpt %d\n", ijpt);exit(1);
    }   
    
    PlotTheo(kFALSE, grstat, grsys, grtheo, itype, ijpt);

    //==========================
    DrawParticle(ipad, nx*2, itype, 1, IDy()==7? (ipad%2==0?0.3:0.75):-999, IDy()==7?0.8:-999);
    //DrawJet(ipad,nx, ijpt, IDy()==7? 0.7:0.3, IDy()==7? (ipad>=nx?0.1:0.3): ipad<nx?0.3:0.10);
    DrawJet(ipad,nx, ijpt, IDy()==7? 0.65:0.3, IDy()==7? (ipad>=nx?0.1:(fgrid==17?0.32:0.3)): ipad<nx?(fgrid==16?0.32:0.3):0.10);
    Draw2Pre(ipad,nx, IDy()==7?0.05:-999, IDy()==7? 0.55: 0.02, IDy()==6?-0.17:-999);
    Draw2Ana(ipad,0, IDy()==7? 0.05 : -999, IDy()==7? 0.0 : (IDx()==2?-0.5:-0.45));
  }

  c3->cd();
  c3->Draw();
  c3->Print(Form("%s_ID%d.eps",c3->GetName(), fgrid));
  c3->Print(Form("%s_ID%d.root",c3->GetName(), fgrid));
  }
  else {
    printf("bad fmode %d\n",fmode); exit(1);
  }
  //_____________________________________________________________________
  //_____________________________________________________________________

  TFile *fout=new TFile(Form("grs%d.root",fgrid),"recreate");
  flist->Write();
  fout->Save();
  fout->Close();

  printf("donedone");
  return 0;
}
